Mott transitions in the half-filled SU(2iVf ) symmetric Hubbard model 
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The Hubbard model with large orbital degeneracy has recently gained relevance in the context of 
ultracold earth alkali-like atoms. We compute its static properties in the SU(2Af) symmetric limit 
for up to M = 8 bands at half filling within dynamical mean-field theory, using the numerically 
exact multigrid Hirsch-Fye quantum Monte Carlo approach. Based on these unbiased data, we 
establish scaling laws which predict the phase boundaries of the paramagnetic Mott metal-insulator 
transition at arbitrary orbital degeneracy M with high accuracy. 
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I. INTRODUCTION 

The interaction-induced Mott transition between a 
metal and a paramagnetic insulator is central to the field 
of strongly correlated electron systems.^ Much insight 
into this phenomenon has been gained in numerical stud- 
ies of the single-band Hubbard model within dynamical 
mean-field theory (DMFT).^ In particular, the phase dia- 
gram and the behavior of characteristic observables (such 
as the effective mass) have been established with high 
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despite the lack of analytic solutions. 



precision' 

While the single-band assumption is rather crude in 
correlated solids (see below) , it can be quite accurate for 
two-flavor mixtures of ultracold fermions on optical lat- 
tices. Since, in addition, the effective interaction between 
neutral ultracold alkali atoms (in their electronic ground 
state) is very short-ranged, such systems appear as nearly 
perfect finite-size realizations of the single-band Hubbard 
model, with the prospect of addressing some of the open 
questions (e.g. regarding high- Tc superconductivity) via 
the tunable quantum simulation of the underlying Hamil- 
tonian. An important step in this direction was the re- 
cent experimental verification of the Mott transition in 
cold-atom systems, ^^'^'^ for which accurate quantitative 
predictions based on DMFT were essential. 

The low-energy electronic properties of correlated 
solids arc usually determined by d orbitals, which are 
five-fold degenerate (per spin) in the atomic limit. This 
degeneracy is partially lifted by the crystal field, result- 
ing, e.g., in a three- fold degenerate t2g band and a two- 
fold degenerate eg band for cubic symmetry. Each of 
these bands is characterized by a local potential plus 
various Hund rule couplings (which can also couple in- 
equivalent bands). Thus, the multi-orbital case is not 
only richer physically-'^^^^^ (including the possibility of 
orbital-selective Mott transitions^^"^^), but is complex 
already by the number of parameters. In addition, ob- 
taining accurate numerical results rapidly becomes costly 
and challenging with increasing number M of orbitals. 
In fact, some methods such as the numerical renormal- 
ization group (NRG) become impractical beyond M = 
2 orbitals. As a consequence, few properties of the 
multi-orbital Hubbard model can be considered well- 
established with high precision, even at the DMFT level. 



However, there is a unique generalization of the single- 
band Hubbard model to arbitrary degeneracy which 
avoids the introduction of additional parameters: In the 
SU(2M) symmetric Hubbard model, all spins and or- 
bitals are equivalent, i.e., share the same local poten- 
tial, the same hopping matrix elements, and they are 
coupled by the same local interaction. In other words, 
the phase space of this particular multi-orbital model 
is identical to the single-band case. Moreover, interest- 
ing analytic insights have been obtained in the limit of 
large band multiplicity M — > c», including an exact ex- 
pression for the critical interaction of the ground-state 
metal-insulator transition (at half band filling) as well 
as scaling arguments for the finite-tcmpcrature critical 
end point. Thus, the sequence of models obtained by 
varying M connects two well-established - and somewhat 
special - limits (M = 1, M = oo), while the interme- 
diate regime M = 2,3,... shares many characteristics 
with generic multi-orbital models, including numerical 
difficulties. Indeed, the SU(2A/) Hubbard model has, 
so far, been explored in this regime only using approxi- 
mate methods, namely the dynamical slave-rotor approx- 
imation (DSR),^^ the projective self-consistent method 
(PSCM),^^'^^ and the self-energy functional approxima- 
tion (SFA) with one bath site per orbital. A fully 
controlled treatment is clearly desirable on fundamental 
grounds and as a solid starting point for generic multi- 
orbital physics. 

Quite recently, the SU(Af) Hubbard model (with to- 
tal degeneracy N > 2) has also become of direct physi- 
cal relevance, namely in the ultracold atom context: In 
rare-earth atoms, a large number of internal states can 
be addressed, which are essentially decoupled from the 
valence electrons. Consequently, all atoms in the elec- 
tronic ground state experience the same optical potential 
and have the same pairwise interactions;^*'^^ a mixture 
with N internal states on an optical lattice can, therefore, 
realize the SU(A^) symmetric Hubbard model. A Mott 
insulating state has already been observed in a SU(6) 
symmetric system of fermionic ytterbium atoms (^^'^Yb) 
on a cubic optical lattice, opening the door to detailed 
experimental investigations of Mott metal-insulator tran- 
sitions in SU(iV) symmetric Hubbard models (with N > 
2). This breakthrough has sparked theoretical inter- 
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est in both SU(A'^) Hubbard^ ^"'^^ and Heisenberg^^ sys- 
tems, with initial studies being limited to one spatial 
dimcnsion^'^^'^^ and to a slave-particle method, respec- 
tively. 

In this work, we construct the phase boundaries of the 
Mott transition at half filling and for up to M = 8 bands, 
based on numerically exact multigrid Hirsch-Fye quan- 
tum Monte Carlo'^^''^^ estimates of characteristic obscrv- 
ables. We also derive scaling laws which predict the phase 
boundaries for arbitrary orbital degeneracy 1 < A/ < oo 
with high accuracy. 

In Sec. II, we establish our notation and relate the 
SU(A^) symmetric Hubbard model to generic multi-band 
models. We also introduce the DMFT in the present con- 
text, discuss our choice of lattice type and energy scales, 
and characterize our DMFT impurity solver. In Sec. Ill, 
we specify the procedure to determine the phase bound- 
aries, briefly summarize literature data for the SU(2A'/) 
symmetric Hubbard model, and present numerically ex- 
act results for M = 2, 4, and 8. Based on these data, 
we deduce in Sec. IV the universal scaling of the critical 
parameters with spin and orbital degeneracy and estab- 
lish the collapse of finite-M data onto an universal phase 
diagram. 



II. MODEL AND METHODS 

The general Hubbard model for M equivalent elec- 
tronic orbitals (e.g., M ~ 3 t2g orbitals) with nearest- 
neighbor hopping t and SU(2) invariant Hund's rule cou- 
pling J has the form: 
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Here, the first line can be viewed as M versions of the 
regular 1-band Hubbard model with (intraorbital) on-site 
Hubbard interaction U; (ij) denotes pairs of nearest- 
neighbor sites i and j, a E {tii} the spin. The cou- 
pling between these orbitals is provided, in general, by 
the interorbital density-density interaction U' and the 
Hund's rule coupling J; the latter contains both an 
Ising-type contribution, coupling to the spin densities 
"»m^ = 4ma'^iraa (sccoud line), as wcll as spin-flip and 
pair- hopping terms (third line). In the limit J — > 0, the 
relation U = U' — 2J implies that the inter- and intraor- 
bital Hubbard interaction become equal: U' — )• U. Thus, 
at J = 0, spin and orbitals are fully equivalent and one 
arrives at the SU(A^) symmetric Hubbard model with 



N = 2M even: 

N 

a=l (ij) 



t^E (2) 

a<Q' 



where a is a combined spin-orbital index. This is pre- 
cisely the situation which has been realized, within the 
single-band approximation and up to the confining poten- 
tial, using rare-earth (earth alkali-like) ytterbium atoms 
{^T^Yh) on a simple cubic optical lattice. Note that 
exchange terms as appearing in (1) at J 7^ require 
a unique classification of the internal degree of freedom 
a S {1, A''} in terms of the "spin" variable a € {tii} 
and cannot arise in the SU(N) symmetric case, where all 
values of a are fully equivalent. 

Paramagnetic Mott metal-insulator transitions (MITs) 
can be expected for this model at all integer fillings 
n = E„("»«> e {l,2,...,Ar- 1} (while n = 0, n = iV 
correspond to band insulators). For N = 2M being even 
(as always in the electronic context) this includes the 
case of half filling n = N/2, where one then expects the 
largest critical interaction (compared to pairwisc equiv- 
alent MITs at fillings n = N/2 ± 1, n N/2 ±2, ... ).^'^ 
In the cases of odd > 3, not to be considered in 
this paper, the Mott plateaus at fillings n = N/2 ±1/2 
are separated by a metallic phase with unique "semi- 
compressible" properties. '^^ 

The DMFT reduces the lattice problem (2) to a 
single-impurity problem,^ with the same local SU(A^) in- 
variant interaction terms, which has to be solved self- 
consistently.'^'' For homogeneous phases, the lattice prop- 
erties enter only via the corresponding tight-binding den- 
sity of states p{e). In line with previous studies, we 
choose the semi-elliptic form associated with the Bethe 
lattice'^® and set the energy scale as ty^ = 1 (for coor- 
dination number Z), which implies unit variance of the 
density of states: dee^p{e) ~ 1. Our numerical re- 
sults can be translated, e.g., to the cubic lattice (in units 
of the hopping t) by multiplying interactions, energies, 
and temperatures by \/Z ~ -\/6 « 2.45.'^^ 

As the interaction couples only to spin-orbital densities 
[i.e., spin flip and pair hopping terms, as arising in the 
general model (1) for J 7^ 0, are absent], DMFT solutions 
can be obtained using quantum Monte Carlo (QMC) im- 
purity solvers without any sign problem for arbitrary 
density. The Hirsch-Fye algorithm'*"''*^ discretizes the 
imaginary-time path integral expression for the Green 
function into A time slices of uniform width At = f3/A, 
where /3 = 1/T (for fce = 1); a Hubbard-Stratonovich 
transformation replaces the electron-electron interaction 
(for each pair a < a') at each time step by a binary 
auxiliary fleld which is sampled using standard Markov 
Monte Carlo techniques. In this work, we use a multi- 
grid implementation'^^''^^ and, thereby, demonstrate that 
its inherent elimination of Trotter errors from the Green 
function and from observables works reliably and accu- 
rately even for a large number M of bands and M (2Af— 1) 
Hubbard-Stratonovich fields. Consequently, our results 
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are free of significant systematic bias, i.e., exact within 
statistical error bars. These statistical errors are reduced, 
compared to a generic M-band model, by employing the 
SU(2Af) symmetry, i.e., by averaging Green functions 
and related observables over all 2M values of the internal 
degree of freedom a and the double (or pair) occupancy 
over all M{2M ~ 1) pairs a < a' . 



III. DETERMINATION OF MIT PHASE 
BOUNDARIES 

It the noninteracting limit — ^ 0, the Hamiltonian (2) 
reduces to the corresponding tight-binding model; due to 
the degeneracy, the system is then metallic at all densi- 
ties < n < A^. In contrast, the energy levels become 
discrete in the atomic limit i ^ 0; at integer filling, the 
system is then an insulator. The question of how the 
evolution between these two limits takes place, e.g. as a 
function of varying U at constant t, has been a matter of 
debate for a long time.^''^'^^'*'^ It is now well-established 
that the (paramagnetic) metallic and paramagnetic in- 
sulating phases are separated, at low temperatures and 
within DMFT, by a sharp transition line in the single- 
band case (i.e., for TV = 2M = 2). 

This transition is of first order at temperatures < 
T < T* (thick blue line in the inset of Fig. 1), evolv- 
ing to second order both at the critical end point (T*, 
U*) and in the hmit T ^ (and U — > Uc2); here and 
in the following, we use the notation Uc2 = Uc2{T ~ 0) 
and Uci — Uci{T = 0) for ground-state values. Due to its 
mean-field character, the DMFT self-consistency equa- 
tions do not directly yield the critical line Uc{T); instead, 
one finds, at T < T*, coexistence of metallic and insu- 
lating solutions in the range Ud (T) < U < Uc2 [T] (indi- 
cated by circles in the inset of Fig. 1). The determination 
of Uc{T) within these boundaries requires a comparison 
of free energies, which arc not directly accessible in QMC 
based approaches (but can be obtained via integration of 
thermodynamic relations"'^). 

As discussed in Sec. I, the situation is expected to be 
quite similar in the multi-band case M > 1. Specifically, 
DMFT should yield a coexistence region of metallic and 
insulating solutions at arbitrary M, including the limit 
M — >■ oo. In this limit, Uc2 was shown^'' to approach 
4|_Eo|j where Eq is the noninteracting ground-state (ki- 
netic) energy; for the Bethe lattice, Eq = — 8M/(37r) « 
— 0.85M. However, as the single-band case deviates from 
this asymptotic value by nearly a factor of two, numerics 
at 1 < A/ < cx) arc needed in order to derive quanti- 
tative predictions from this analytic result. Linearized 
DMFT^^ is not sufficient in this respect; its prediction 
Uc2 = 4M -I- 2 (for arbitrary M) is consistent with the 
exact asymptotic result only regarding the power (in M), 
but the prefactor (4 instead of 32/(37r) sa 3.4) is ob- 
viously incorrect. Remarkably, the critical interaction 
at finite temperatures scales differently: U* oc M-'^/^, 
as was argued convincingly analytically.^^ In this case. 
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FIG. 1: (Color online) Previous results for the Mott metal- 
insulator transition in the SU(2M) symmetric Hubbard model 
within dynamical mean-field theory: coexistence phase dia- 
grams for band degeneracy M = 1, 2, and 4 obtained us- 
ing the dynamical slave-rotor approximation (DSR)^^ and the 
self-energy functional approximation (SFA),^^ respectively, in 
comparison with numerically exact quantum Monte Carlo 
(QMC) data for M = 1. Inset: magnified view for M = 1. 



the analytic considerations do not even yield a prefac- 
tor; thus, quantitative predictions regarding U* are com- 
pletely dependent on accurate numerical results for suf- 
ficiently large values of M . 



A. Previous results for M > 2 

So far, the complete DMFT coexistence regions have 
been computed at M > 1 only using the dynamical 
slave-rotor formalism^^ and using the self-energy func- 
tional approach. The former is an approximate impu- 
rity solver which contains a free parameter and had been 
tested quantitatively only for M = 1. Consequently, the 
accuracy of its results at M > 1 (and even in the limit 
M — > oo) is, a priori, completely unclear. In contrast, the 
SFA^^ is based [in the variant used in Ref. 27, known as 
the dynamical impurity approximation (DIA)] on a dis- 
cretization of the DMFT dynamical bath; it reduces to 
the DMFT in the limit of an infinite number of bath sites, 
A^h — > oo. So this method is numerically exact (within 
DMFT); however, an unknown bias remains for a finite 
value of N\j, in particular for the "two-site SFA" with a 
single bath site (per interacting orbital), as employed in 
Ref. 27. It is, a priori, unclear, how this bias evolves with 
M (at fixed Nb/M). 

As shown in Fig. 1, the DSR (dotted lines) and the SFA 
(dashed lines) both yield coexistence regions for M = 2 
and A/ = 4 which have shapes similar to those in the 
single-orbital case M = 1. At M = 2, even the critical 
interactions are in good mutual agreement with a value 
U* ~ 6.3; however, this agreement seems coincidental, as 
the DSR estimate of U* is significantly below (above) the 
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SFA estimate at M = 1 (M = 4). In general, the DSR 
appears to yield much larger coexistence regions than the 
SFA. As the DSR is an uncontrolled and comparatively 
cheap approximation, one might be tempted to put more 
trust in the SFA results. However, both the DSR and 
the SFA deviate very significantly from the exact QMC 
results previously established for M = 1: as seen in the 
inset of Fig. 1, the SFA underestimates the critical tem- 
perature T* by about a factor of 2 and the area of the 
coexistence region by even more, while the DSR overes- 
timates the latter by nearly the same factor. Given these 
discrepancies, it is clear that the (previously published) 
data shown in Fig. 1 are not sufficient for verifying the 
scaling laws discussed above and for determining their 
prefactors and corrections at finite M . 



B. Insights from the single-band case {M = 1) 

In order to achieve this goal, we will, in the remain- 
der of this section, determine unbiased coexistence phase 
boundaries at M = 2 and M = 4 and determine T* and 
U* at M ~ 8, based on exact QMC data. For complete- 
ness and for better illustration of the asymptotic behavior 
of the relevant observables in the limit T — > 0, we will 
first discuss results for the single-band case (Af = 1), 
depicted in Fig. 2. 

The quasiparticle weight Z = m/m* quantifies the 
rcnormalization of the quasiparticlcs in a Fermi liquid 
by interactions and is closely associated with the inverse 
linear specific heat: Z{U,T = 0) 7(0)/7(f7) [with en- 
ergy E{U, T) = E{U, 0) + j{U) + 0{T^)]. It can be 
expressed (exactly) in terms of the self-energy S(a;) as 



(a) 0.: 



Z-^ = 1 - dRe^{uj)/duj\ 



(3) 



Fig. 2(a) shows corresponding discrete QMC estimates at 
finite temperatures, based on the value of the self-energy 
at the first Matsubara frequency iuji = inT: 



1 -ImS;(i7rr)/(7rr). 



(4) 



Clearly, the data set for each temperature (denoted by 
symbols) is split into two branches: one metallic branch 
with moderately high values {Z > 0.04) which extends 
down to U = (shown only for U > 4.5) and an insu- 
lating branch where Z « for the lower temperatures 
(e.g. T = 0.02, denoted by upward triangles). Only at 
the highest temperatures shown (T = 0.04 and T = 0.05) 
do the estimates in the insulating phase reach values up 
to about 0.01, which is mainly an artifact of the discrete 
approximation of Z, Eq. (4); in particular, these values 
of Z have no relation to the specific heat (which is expo- 
nentially small in this range). 

In contrast, the double occupancy D = (nifriii) (i.e., 
the probability of a site being occupied by two fermions 
simultaneously), which is depicted in Fig. 2(b) and re- 
lated to the interaction energy by i^int = UD, has very 
significant values in both phases. D is independent of 
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FIG. 2: (Color online) Numerically exact DMFT results (sym- 
bols), obtained using multigrid HF-QMC, in the vicinity of 
the Mott metal-insulator transition for M = \: (a) quasipar- 
ticle weight (b) double occupancy (= pair occupancy) D, 
and (c) energy i? as a function of the on-site interaction U . 
The extrapolations to the ground state (black solid lines) for 
D and E include perturbative informations and thermody- 
namic consistency;^^ other lines are guides to the eye only. 



temperature at the scale of the figure in the insulating 
phase, within its temperature-dependent range of stabil- 
ity, i.e., for U > C/ci(r). This behavior is expected in a 
gapped phase, where thermal excitations are suppressed 
exponentially. It allows high-precision estimates of the 
ground-state function Dins{U) from QMC, in particu- 
lar by extrapolation of high-order coefficients of strong- 
coupling perturbation theory.^^ On the metallic side, D 
depends strongly on temperature, especially in the range 
T < T* « 0.055.4'' function of U, the shapes of 

the curves look remarkably similar for D and Z; for 
both observables, the (negative) curvature becomes much 
stronger near the boundaries of the metallic phase, i.e., 
at [/ < Uc2iT). 

In comparison, the results for the energy E = (H) 
look nearly linear in Fig. 2(c) as function of U in the 
same parameter ranges (which implies that the kinetic 
energy, not shown, has a positive curvature which nearly 
cancels that of D); the values also approach those of 
the insulating solution (again with invisible tempera- 
ture dependence) much more closely. As the relation 
D = dF/dU for free energy F = E — TS reduces to 
Dmct{U) = dEmct{U)/dU in the ground state, D{U,T) 
and £'([/, T) are not independent at low T; this connec- 
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tion as well as the relation between 7 and Z have made it 
possible to determine the ground-state energetics [black 
solid lines in Fig. 2(b) and 2(c)] in the metallic phase as 
wen." 

The crucial point for determining the phase boundaries 
[/ci(r), Uc2{T) via data sets such as depicted in each of 
the panels in Fig. 2 is that all included data points actu- 
ally denote converged solutions, i.e., they correspond to 
fixed points of the DMFT self-consistency cycle, whereas 
no metallic solutions exist at U > Uc2{T) and no insu- 
lating solutions exist at f/ < Uci{T), respectively. Both 
the verification and the exclusion of such fixed points 
are very difficult to achieve reliably, as numerical noise 
(associated with Monte Carlo importance sampling for 
a finite number of sweeps), systematic bias (e.g. result- 
ing from Trotter errors) and critical slowing down (for 
T w T* and U ~U*) can easily lead to false positives or 
negatives. For this reason, it is essential to monitor sev- 
eral observables at the same time, as deviations from the 
expected systematics can help to identify artifacts of in- 
complete convergence or divergence after a (necessarily) 
finite number of DMFT iterations. 

In this manner, we have obtained the phase bound- 
aries shown as circles in the inset of Fig. 1 (building 
upon earlier work^^) with high precision at finite tem- 
peratures T > 0.01. The squares denote complementary 
ground-state results for Ud from extrapolated perturba- 
tion theory^-^ and for Uc2 from ED and NRG.^'* Taken 
together, these results determine fit functions for the co- 
existence region (thin solid lines and blue-shaded region) 
with high precision; we will later test the hypothesis that 
very similar fits might capture the Mott transition at 
M > 1. Note that the (numerically exact) QMC es- 
timate of Uc2{T) and its extrapolation to T — > agree 
well with the corresponding SFA estimate (at T < Tgp^). 
The inset of Fig. 1 also shows a thick line within the co- 
existence region that denotes the DMFT estimate of the 
actual first-order phase transition.*^ 



C. QMC results for M > 2 

The quasiparticle weight, as defined above, remains 
a well-defined and useful concept at arbitrary degener- 
acy and is shown in Fig. 3(a) for M = 2. However, 
as more than two fermions can occupy the same site 
for A/^ > 2, it is advantageous to generalize the concept 
of the double occupancy to that of the pair occupancy 
D = X]q<q' (^ia'^ict') j have retained the symbol "D" 
for this observable, as it obviously reduces to the dou- 
ble occupancy in the single-band case and satisfies the 
relation Ei^t = DU, stated in the previous section, for 
arbitrary degeneracy N (or number of orbitals M). At 
fixed integer band filling n, its minimum value as a func- 
tion of U and T is -Dmin = n{n~l)/2, corresponding to an 
atomic state with exactly n filled orbitals; while this min- 
imum is zero in the single-band case at half filling (n = 1), 
it has the values 1, 6, and 28 in the half-filled case 
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FIG. 3: (Color online) Numerically exact DMFT-I-QMC re- 
sults (symbols) for AI = 2: (a) quasiparticle weight Z, (b) 
pair occupancy D, shifted by its value in the atomic limit, 
(c) energy E, relative to the asymptotics in the atomic limit; 
lines are guides to the eye only. 



{n = M) at M = 2, M = 4, and M = 8, respectively.'''^ 
For better comparison to the previous results, we have, 
therefore, subtracted -Dmin = 1 in Fig. 3(b). The corre- 
sponding interaction energy i^min = U -Dmin = U has also 
been subtracted from the energy in Fig. 3(c); only with 
this adjustment does the slope dE/dU approach zero in 
the limit of strong interaction (in the insulating phase). 

With these adjustments, the data shown in Fig. 3 for 
M = 2 look remarkably similar to the single-band case 
at low temperatures; in addition to data for T < T*, 
i.e., with coexistence, we have also included results for 
T = 0.118, where the DMFT solution is unique for all in- 
teractions, corresponding to a continuous crossover curve 
(solid line). 

By reading off the phase boundaries from these numeri- 
cally exact QMC data we can, for the first time, construct 
the coexistence phase diagram for M = 2 in an unbiased 
way, down to the lowest QMC temperature T = 0.04, 
as denoted by circles in Fig. 4. Also shown is the SFA 
prediction^^ (green dashed lines and green-shaded area) 
as well as the DSR result^^ (dotted lines) . Quite remark- 
ably, the QMC estimates of Uc2 agree perfectly with the 
corresponding SFA prediction at T < Tgp^ (for M = 2), 
even better than in the case M = 1 (cf. inset of Fig. 
1). Consequently, we regard two-site SFA as practically 
exact (only) for Uc2{T) at M > 2 and will not try to 
compete with its estimates for the corresponding ground- 
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FIG. 4: (Color online) Coexistence phase diagram for M — 2 
(within DMFT): exact QMC resuhs (circles) for the bound- 
aries at T > 0.04 are consistent with a fit (solid lines and 
shaded area) obtained by rescaling the exact M = 1 coex- 
istence region. Also shown; SFA result^^ (dashed lines and 
green-shaded area) and DSR prediction^^ (dotted lines). 



state value Uc2- In contrast, the exact QMC results for 
Uci{T) are significantly lower than their SFA counter- 
parts; also the QMC value for T* « 0.12 is significantly 
above the SFA estimate. At the same time, the QMC 
data for M = 2 (circles) are in excellent global agreement 
with a rescaled version (solid blue lines) of the numeri- 
cally exact one-band result (solid blue lines in the inset of 
Fig. 1) determined above; thus, the coexistence regions 
for M = 1 and M = 2 appear to be similar even in the 
strict mathematical sense. 

The DSR yield phase boundaries (dotted lines) of very 
similar shape, but shifted towards larger U and T, with a 
prediction for U* which nearly coincides with that of the 
SFA. Our exact data now reveal that both estimates are 
too large by about 0.3. Still, both approximate methods, 
SFA and DSR, yield more accurate predictions of the 
phase diagram for M = 2 than for M = 1 (which in the 
case of the DSR might be due to a specific parameter 
choice); in particular, the discrepancies with respect to 
the area of the coexistence region (of about 20%) are 
much smaller. 

The QMC results for M = 4, shown in Fig. 5, include 
three temperatures very close to T*: both at T = 0.25 
(squares) and at T = 0.222 (circles) the curves are contin- 
uous, without coexistence, while a clear coexistence is ob- 
served at T = 0.2 (upward triangles). At the same time, 
the maximum derivatives dZ/dU and dD/dU are much 
larger at T = 0.222 than at the neighboring grid points; 
we conclude that T* w 0.22. QMC results, which are 
already much more expensive computationally at M — 4 
than in the single-band case [the cost being roughly pro- 
portional to the number M(2M - 1) = 28 of Hubbard- 
Stratonovich fields] have also been obtained near T* /2. 
Overall, the evolution of all three observables (Z, D, and 
E 3S, a. function of U and T) is consistent with the ex- 
pectations from M = 1 and M = 2 (cf. Fig. 2 and Fig. 3, 
respectively) within symbol sizes, which reflect approxi- 
mate error bars. 

Corresponding phase boundaries are shown as circles 
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FIG. 5: (Color online) Numerically exact DMFT-fQMC re- 
sults (symbols) for M = 4, analogous to Fig. 3. 
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SFA result^^ (dashed lines and green-shaded area) and DSR 
prediction^^ (dotted lines). 



in Fig. 6. These data confirm again, both the accuracy 
of the SFA prediction for Uc2 and the validity of the scal- 
ing assumption for the shape of the coexistence region, 
yielding the blue solid lines in Fig. 6. 

The scaling assumption is also supported by the ob- 
served convergence of the SFA results for Ud (T) towards 
these rescaled one-band results (blue solid lines) in the 
series M = 1 (inset of Fig. 1), M = 2 (Fig. 4), and M = 4 
(Fig. 6); in the last case, the SFA discrepancy in T* has 
already shrunk to about 5% and that in the area of the 
coexistence region to about 10%. While the DSR yields 
an essentially correct value of T* at M = 4, the whole 



7 



Q 0.1 




T= 0.333 -B- 
T= 0.308 -0- 
T= 0.286 - A- 
J I 



FIG. 7: (Color online) Numerically exact DMFT+QMC re- 
sults (symbols) for AI = 8, analogous to Fig. 3 and Fig. 5. 
Increased symbol sizes reflect larger error bars. 



DSR coexistence region appears shifted towards larger 
interactions (relative to the unbiased QMC results, de- 
noted by circles and solid lines) by an offset of roughly 
1/8 of its true width, similarly to the case M = 2. We 
conclude that the DSR, in contrast to the SFA (with one 
bath site per orbital), does not become more accurate at 
large AI. 

Let us, finally, turn to the case of M = 8 orbitals (i.e., 
a total degeneracy of 16 in the spin-|-orbital space), which 
has never been considered in the literature before. Due 
to the extreme computation cost (increased by a factor 
of 120 relative to the 1-band case), and since we have 
already established the universal shape of the phase dia- 
gram, we focus on temperatures in the immediate vicinity 
of the critical point. The QMC results depicted in Fig. 7 
show an increased scatter, indicating larger relative error 
bars as represented by the increased symbol sizes. 

Still, they allow us to locate the critical point at T* w 
0.33, U* w 10.9. Together with the value Uc2 ~ 29.6 
read off from Fig. 9, these parameters also determine 
Uci « 12.8 and (using the known shape) the full coex- 
istence phase for M ~ 8, shown in Fig. 8. Due to the 
exponential scaling of exact diagonalization with the to- 
tal number of orbitals (interacting and bath), an SFA 
solution analogous to those shown for M = 1,2,4 would 
be prohibitively expensive at M = 8. 



IV. SCALING OF CRITICAL PARAMETERS 
WITH SPIN AND ORBITAL DEGENERACY 

As seen in the preceding section, the critical parame- 
ters T* , U* (of the finite-temperature critical end point) 
and Uc (of the ground-state Mott transition) all increase 
significantly with increasing degeneracy [i.e., with the 
number M of orbitals, corresponding to = 2M fold 
degeneracy in the spin-|-orbital space]. In order to study 
the dependencies in detail, the left column of Fig. 9 shows 
the estimates of these parameters as a function of the in- 
verse number of orbitals, M~^. At the scale of Fig. 9 (b), 
all finite-temperature methods give quite similar results 
for U* , with slight deviations for DSR (triangles and dot- 
ted line) at large degeneracy. Deviations become much 
more apparent for T*, shown in Fig. 9 (c), with SFA data 
(diamonds and dashed line) having a nearly constant neg- 
ative offset relative to the exact QMC data (circles and 
solid line). Regarding Uc2, we see in Fig. 9 (a) that DSR 
is far above the accurate SFA data at larger M; we have 
also included the L-DMFT prediction Uc2 = 4M + 2. Ob- 
viously, all of the observables increase strongly towards 
smaller 1/M (i.e., towards larger degeneracy); alas, it is 
hard to distinguish exponents at this level. 

The scaled data shown in the right hand column of 
Fig. 9 demonstrate convincingly, however, that Uc2 in- 
deed scales as M while U* scales as M^/^; in addition, 
we establish that also T* scales as M^/^. In Fig. 9 (d) we 
have, specifically, divided Uc2 by 2M -|- 1 (instead of M) 
in order to convert the L-DMFT prediction to a constant 
(with value 2). With this particular scaling ansatz, also 
the SFA data fall on a straight line, interpolating between 
the numerically exact results for Af = 1 and the analytic 
expression for M — oo (square). Our fit corresponds to 
the scaling law 



U^2 ~ 1-70 (2A/ + 1) (1 + 0.166 A/"^) 



(5) 



In this representation, DSR is even seen to have the 
wrong tendency; this method should be off by more than 
a factor of two for Af — )■ oo . 

The same offset in the argument is also seen, in Fig. 
9 (e), to minimize curvature when rescaling estimates 
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of U* (to U*/V2MT1). Specifically, the exact QMC 
data become nearly flat and can safely be extrapolated 
to 1/M = 0, with the result (given without higher order 
corrections as they are not significant) 



U* w 2.67V27\/ + 1 



(6) 



The SFA data are significantly above the QMC results 
at all finite M. In the extrapolation to 1/M = some 
discrepancy remains; it is not entirely clear whether it is 
significant. 

For rcscafing T*, we have chosen a different offset in 
Fig. 9 (f) which ensures, again, that the results of each 
method fall on a nearly straight line, with perfect con- 
vergence of the SFA data to the exact QMC results. We 
conclude that T* is well represented by the expression 



T* w 0.090 V2M - 1 (l - 0.41 M' 



(7) 



Note that both corrections to the asymptotic scaling 
T* (X M^/^, arising from the shift in the argument and 
associated with the explicit 1 /M term, work in the same 
direction: in the physical range of Af , the critical temper- 
ature increases much faster with the degeneracy than one 
would expect from the scaling law. For example, going 
from SU(2) to the SU(6) Hubbard model, recently real- 
ized with ultracold rare-earth atoms, increases T* by a 
factor of 3.3, much larger than the factor sa 1.73 sug- 
gested by the large-M asymptotics. This extra enhance- 
ment is certainly beneficial for accessing Mott physics in 
cold atom experiments. 
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FIG. 10: (Color online) Scaling phase diagram: upon rescal- 
ing with the parameters T* , U* of the second-order finite- 
temperature critical end point and with the critical interac- 
tion Uc2 of the second-order ground-state Mott transition, the 
exact QMC data (symbols) collapse onto the scaling curves 
(solid lines). SFA results (broken lines and shaded regions) 
deviate, but converge towards this scaling form for large M. 
Inset: DSR results^* (broken lines) deviate from the scaling 
form (solid lines) at all M; in contrast, a high-precision SFA 
calculation''* (with 5 bath sites: circles) nearly recovers the 
reference result already at M = 1. 



Let us stress again, that all numerical results corre- 
spond to a semi-elliptic density of states of unit variance 
(and full bandwidth 4); they can be converted to the 
cubic lattice, in units of the hopping t, by multiplication 
with y/6 sa 2.45 (or to the square lattice by multiplication 
with Vi = 2).39 

The expressions (5) - (7) fully determine the coexis- 
tence phase diagram at any orbital degeneracy M when 
combined with the scaling phase diagram Fig. 10 in 
which, by construction, the finite-temperature critical 
point has the coordinates (0,1) while the ground-state 
critical point has the coordinates (1,0) for any value of 
M. Its main panel shows that the QMC data (symbols) 
for the phase boundaries indeed collapse onto a universal 
phase diagram (black solid lines and gray shaded back- 
ground) in this representation, while the SFA data (using 
one bath site per interacting orbital) approach it only at 
large M. As seen in the inset of Fig. 10, the inclusion of 
a larger number of bath sites (instead of one per orbital) 
vastly improves the accuracy of the SFA also at M = 1 
(filled circles),*® beyond the level of the multi-band case 
with the same total number of sites. The inset further 
shows that the DSR data seem to converge after rescaling 
(with a ncar-coUapsc between the results for M — 2 and 
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M = 4), but to an incorrect limit. 

From the universal phase diagram, one can also read 
off that the insulating state is meta-stable (within DMFT 
and in a paramagnetic phase) at zero temperature down 
to 

0.9(7* + 0.1 C/e2; (8) 

which is easily expressed cxplicitcly in terms of M by 
using Eqs. (5) and (6). 

V. CONCLUSION AND OUTLOOK 

In conclusion, we have studied the Mott metal- 
insulator transition of the SU(2A/) symmetric Hubbard 
model by solving the paramagnetic DMFT equations nu- 
merically exactly. Our results confirm the predicted^^ 
asymptotics of the ground-state critical interaction Uc2 oc 
M for M — >■ oo and determine the (previously un- 
known) subleading corrections. They also confirm the 
predicted^^ exponent (of 1/2) of the dependence of the 
finite-temperature critical interaction U* on M and yield 
the missing prefactor plus subleading terms; in addition, 
they establish the relation T* oc M^/^ also for the asso- 
ciated critical temperature. Despite the different scaling 
of the end points of the first-order phase-transition line 
with M, the shape of the coexistence region is found to 
be universal to an astonishing degree. This universality 
could only be revealed by a method (niultigrid HF-QMC) 
that is numerically exact at arbitrary M; in earlier SFA 
and DSR studies, it was obscured by systematic errors. 

Our results yield precise predictions for the Mott tran- 
sition at arbitrary values of M, to be tested in cold-atom 



experiments. Due to the enhanced critical temperatures, 
the multi-flavor case might make Mott physics more ac- 
cessible than in the single-band (i.e., two- flavor) case, 
in which the Mott signatures^^'^'^ seen so far correspond 
to crossovers, not true phase transitions. On the other 
hand, the experimental two-flavor studies profited from 
the fact that the MIT extends, as a crossover, far above 
T* with relatively little variation in U; thus, it is possible, 
e.g., to obtain good estimates of U* from measurements 
at T > T*. This is stiU true in the SU(3) case.^e At large 
M, however, the relative variation of U along the MIT line 
increases significantly, from {Uc2 — U*)/Uc2 ~ 0.2 at M = 
1 to, e.g., {Uc2 - U*)/Uc2 ~ 0.6 at M = 8.^9 One may 
suspect that the relative variation of U in the crossover 
region is similarly enhanced at large degeneracy, which 
implies that a closer approach of T* would be required 
in order to determine U* . Such low-temperature experi- 
ments might also explore ordering phenomena, which are 
a fascinating topic of their own and beyond the scope of 
this paper. 

More generally, our results provide high-precision 
numerical benchmarks for evaluating DMFT impurity 
solvers in the challenging regime of moderate to high 
orbital degeneracy; they could also be used for assess- 
ing the relative importance of nonlocal correlations at 
higher band degeneracy, e.g. by comparison with high- 
temperature expansions'*^ or with direct exact calcula- 
tions once they become available. 
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